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Abstract 

A geometric derivation of nonholonomic integrators is developed. It is based in the classical 
technique of generating functions adapted to the special features of nonholonomic systems. The 
theoretical methodology and the integrators obtained are different from the obtained in m- 
In the case of mechanical systems with linear constraints a family of geometric integrators 
preserving the nonholonomic constraints is given. 

AMS classification scheme numbers: 37J60, 58F05, 37M15 


1 Introduction 

1.1 Introduction to nonholonomic mechanics 

The theory of systems with nonholonomic constrains goes back to the XIX century. D’Alembert’s 
or Lagrange-D’Alembert’s principle of virtual work and Gauss principle of least constraint can be 
considered to be the first solutions to the analysis of systems with constraints, holonomic or not. 
After a period of decay, recently many authors show a new interest in that theory and also in its 
relation to the new developments in control theory, subriemannian geometry, robotics, etc (see, for 
instance, 0H). The main characteristic of this period was that Geometry was used in a systematic 
way (see L.D. Fadeev and A.M. Vershik ESI as an advanced and fundamental reference and, also, 

El IH IZl CHI DEI ESI HH HZ1 EHl EH ) 

As is well known, in most problems of particle mechanics, the motion of the particles is constrained 
in some way; this is the term used to denote the condition that some motions or configurations are 
not allowed. First, we will start with a configuration space Q, which is a n-dimensional differentiable 
manifold, with local coordinates q l . General two-side or equality constraints are functions of the 
form <t> a {q l ,c?) = 0,1 < a < m, depending, in general, on configuration coordinates and their 
velocities. The various kinds of constraints we are concerned with will roughly come in two types: 
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holonomic and nonholonomic, depending whether the constraint is derived from a constraint in the 
configuration space or not. Therefore, the dimension of the space of configurations is reduced by 
holonomic constraints but not by nonholonomic constraints. Thus, holonomic constraints permit 
a reduction in the number of coordinates of the configuration space needed to formulate a given 
problem (see .11]). 

We will restrict ourselves to the case of nonholonomic constraints, since the case of holonomic 
constraints, and, in particular, the construction of holonomic integrators, is well established in the 
existing literature. Geometrically, nonholonomic constraints are globally described by a submani¬ 
fold M of the velocity phase space TQ, the tangent bundle of the configuration space Q. In case 
M is a vector subbundle of TQ, we are dealing with linear constraints. We will usually refer to 
M as D and, in such case, the constraints are alternatively defined by a distribution D on the 
configuration space Q. If this distribution is integrable, we are precisely in the case of holonomic 
constraints. In case M is an affine subbundle modeled on a vector bundle D, we are in the case of 
affine constraints. In the sequel, we will denote by D the constraint submanifold on the velocity 
phase space, no matter if they are determined by linear or nonlinear constraints. 

Given the constraints, we need to specify the dynamical evolution of the system. The central 
concepts permitting the extension of mechanics from the Newtonian point of view to the Lagrangian 
one are the notions of virtual displacements and virtual work; these concepts were formulated 
in the developments of mechanics, in their application to statics. In nonholonomic dynamics, 
the procedure is given by Lagrange-D’Alembert’s principle. We usually consider nonholonomic 
constraints of linear type, which are the constraints that we will regard as natural in a mechanical 
sense (although the extension for general nonholonomic constraint will be straightforward). We now 
come to the description of the constraint forces; for constraints of that type, Lagrange-D’Alembert’s 
principle allows us to determine the set of possible values of the constraint forces only from the 
set of admissible kinematic states, that is, from the constraint manifold D determined by the 
vanishing of the nonholonomic constraints. Therefore, assuming that the dynamical properties of 
the system are mathematically described by a configuration space Q, by a Lagrangian function 
L and by a distribution determining the linear constraints D , the equations of motion, following 
Lagrange-D’Alembert’s principle, are 



(1) 


where 5q l denotes the virtual displacements verifying 


rtW = 0 


( 2 ) 


and D° = span {/i“ = (for the sake of simplicity, we will assume that the system is not 

subject to non-conservative forces). By using the Lagrange multiplier rule we obtain that 




The term on the right represents the constraint force or reaction force induced by the constraints. 
The functions A a are Lagrange multipliers to be determined in order to obtain a set of second order 
differential equations. These Lagrangian multipliers are computed using the constraint equations. 
An interesting remark, that will be used in the sequel, is that whenever the Lagrange multipli¬ 
ers X a = A a (q l ,q l ) have been determined, then the system of equations © can be considered a 
Lagrangian system subject to external conservative forces given by the right-hand side term, tak¬ 
ing, obviously, an initial condition on the constraint submanifold D. Automatically, the choice of 
the Lagrange multipliers A a implies that the solution integral curves also verifies the constraint 
equations. 

1.2 Introduction to Geometric Integration and Discrete Mechanics 

Standard methods for simulating the motion of a dynamical system, generically called numerical 
integrators, usually take an initial condition and move it in the direction specified by the equation 
of motion or an appropriate discretization. But these standard methods ignore all the geometric 
features of many dynamical systems, as for instance, for Hamiltonian systems we have preservation 
of the symplectic form, energy (in the autonomous case) and symmetries, if any. However, new 
methods have been recently developed, called geometric integrators, which are concerned with some 
of the extra features of geometric nature of the dynamical systems. Usually, these integrators, in 
simulations, can run for long times with lower spurious effects (for instance, bad energy behavior 
for conservative systems) than the traditional ones. As is well known, the typical test example 
is the simulation of the solar system. Therefore, there is presently a great interest in geometric 
integration of differential equations as, for instance, symplectic integrators of Hamiltonian systems 

HU 371 . 

Discrete variational integrators appear as a special kind of geometric integrators. These integrators 
have their roots in the optimal control literature in the 1960’s and 1970’s (Jordan and Polack d, 
Cadzow J5], Maeda jSU S3) and in 1980’s by Lee mad, Veselov muni- In these papers, there 
appear the discrete action sum, discrete Euler-Lagrange equations, discrete Noether theorem... 
Although this kind of symplectic integrators have been considered for conservative systems 
13511121 [50115T| . it has been recently shown how discrete variational mechanics can include forced 
or dissipative systems HD EE], holonomic constraints HU mi, time-dependent systems EM HU, 
frictional contact m and nonholonomic constraints (see pi mi)- Moreover, it has been also 
discussed reduction theory SHIM ED, extension to field theories PESM and quantum mechanics 
M- All these integrators have demonstrated exceptionally good longtime behavior and the research 
of this topic is interesting for numerical and geometric considerations. 

Now, we will describe the discrete variational calculus, following the approach in m (see also 
mm)- a discrete Lagrangian is a map : Q x Q —* M (this discrete Lagrangian may be 
considered as an approximation of the continuous Lagrangian L : TQ —> M). Define the action sum 
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Sd '■ Q N+1 —> ffi corresponding to the Lagrangian L d by 

N 

d ^ ' ^d^Qk—li Qk) j 
k =1 

where qk E Q for 0 < k < N. For any covector a E TA X2 ^(Q X Q), we have a decomposition 
a = ai + «2 where E T*.Q. Therefore, 

dLd{qo, Qi) = DiLd(qo, Qi) + D 2 L d (q 0 , qi) . 

The discrete variational principle or Cadzow’s principle states that the solutions of the discrete 
system determined by L d must extremize the action sum given fixed points qo and qysr■ Extremizing 
Sd over qk, 1 < k < N — 1, we obtain the following system of difference equations 

D\L d {qk, Qk+i) + D 2 Ld{qk-\,qk) = 0 . 

These equations are usually called the discrete Euler-Lagrange equations. Under some regularity 
hypothesis (the matrix (Di 2 Ld(qk,qk+i)) is regular) this implicit system of difference equations 
defines a discrete flow T : Q x Q —> Q x Q, by T(qf.-i, qk) = qk+ 1 )- 

The geometrical properties corresponding to this numerical method are obtained defining the dis¬ 
crete Legendre transformation associated to L d by 

FL d : QxQ —f T*Q 

(qo,qi) 1 —^ (qo,-DiLd{qo,qi)) , 

and the 2-form u> d = FL*,ujq, where loq is the canonical symplectic form on T*Q. The discrete 
algorithm determined by T preserves the symplectic form u>d, he., = ujd- Moreover, if the 

discrete Lagrangian is invariant under the diagonal action of a Lie group G, then the discrete 
momentum map J d ■ Q x Q Q* dehned by {Jd{qk,qk+ 1 ),£) = (D 2 L d (q k , qk+i), Cq(%+i)) is 
preserved by the discrete flow. Therefore, these integrators are symplectic-momentum preserving 
integrators. Here, is the fundamental vector field determined by £ E 0. 

Another alternative approach to discrete variational calculus comes from the classical theory of 
generating functions (see, for instance, [I]). Since (T*Q,u>q) is an exact symplectic manifold, where 
loq is the canonical symplectic form of T*Q and ujq = —d0Q, the symplectic flow Fh : T*Q —> T*Q of 
a Hamiltonian vector field Xu is a canonical transformation, and then GraphfiqJ, the graph of Fy,. 
is a Lagrangian submanifold of the symplectic manifold ( T*Q x T*Q , Ll) where Ll = tt 2 ujq — vrj loq. 
Here, we denote by ny : T*Q x T*Q —>• T*Q , i = 1,2 the canonical projections. Therefore, denoting 
0 = tt^Oq — 7t*0q we have that 

i* Fh Q = ~di* Fh O = 0 , 

where ip h : Graph(FX) i—» T*Q x T*Q is the canonical inclusion. Then, at least locally, there exists 
a function S h : Graph(F) l ) —» M such that i Fh & = dS h . Taking ( q l ,pi ) as natural coordinates in 
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Graph(F/j) and (q l ,Pi, q\ p*) the coordinates in T*Q x T*Q, then, locally S h is a function of (q,p) 
coordinates. Hence, along Graph(F/j), we have q 1 = q l (q,p) and p l = p l (q,p) and moreover 

Pi dq l - pidq 1 = dS h (q,p) . 

Assume that in a neighborhood of some point x € Graph(i 7 ) i ), we can change this system of 
coordinates by new independent coordinates ((?\q*) (the local condition is that det(<9q/<9p) 0). 

In such a case, the function S h can be locally expressed as S h = S h (q,p ) = S h (q, q). The function 
S h (q, q) will be called a generating function of the first kind of the canonical transformation T},. 
Moreover, 

dS h 

Pi ~ ’ 

dS h 


P* = 


dq l 


A nice and useful interpretation of the discrete Euler-Lagrange equations is the following theorem 

[23 EH- 

Theorem 1.1 Let the function S Nh be defined by 

N-l 


S Nh (qo,qN) = sh ( < lk,Qk+ 1 ) , 


k =0 


where qk, 1 < k < N — 1, are stationary points of the right-hand side, that is 
0 = D2S h (qk-i,qk) + DiS h (qk,qk+i) , l<k<N—l 


(4) 


then S Nh is a generating function of first class for F^h '■ T*Q —> T*Q, for h sufficiently small and 
where F^h denotes the flow of Xh over time Nh. 

Moreover, if we start with a regular Lagrangian function L : TQ —> M, and Ft : T*Q —>• R is the 
locally associated Hamiltonian, then we also have the following result (for example, see I32j ) 

Proposition 1.2 A generating function of the first kind for Fu is given by 


S (qo,qi)= L(q(t),q(t))dt , 

Jo 

where q(t ) is a solution of the Euler-Lagrange equations such that q(0) = qo and q(h ) = q\. 

The conclusion is that the discrete variational calculus consists in taking an approximation of 
the generating function S h . From this approximation we obtain a new Lagrangian submanifold of 
T*QxT*Q and the relation between subsequent steps is given by Q for the new generating function, 
which are precisely the discrete Euler-Lagrange equations. The symplecticity and preservation of 
momentum are now direct consequences of this description. 
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1.3 Introduction to nonholonomic integrators 

In a recent paper, J. Cortes and S. Martinez m have proposed a construction of nonholonomic 
integrators which is useful for numerical considerations. Their construction is based on the discrete 
Lagrange-D'Alembert's principle. Assuming that the constraints are given by a distribution D , this 
principle states that 

(■ DiL d (q k ,q k+1 ) + D 2 L d (q k _ 1 ,q k )) i Sq l k =0, 1 < i < N - 1 , 

where 5q k E D qk and, in addition ( q k ,q k +i ) E D d . Here D d denotes a discrete constraint space 
D d C Q x Q. This integrator has a good performance and naturally inherits some geometric 
properties of the continuous problem. Observe that the method is based on the discretization of 
the Lagrangian and a coherent discretization of the constraints, and both determine the discrete 
constraint forces. 

Alternatively, we propose a nonholonomic integrator also based on the discretization of the La¬ 
grangian function (in a more precise sense, we discretize the action function) but now we take 
a coherent discretization of the constraint forces and both determine the discrete constraint sub¬ 
manifold. This method gives us, in general, different integrators from those in m- The last 
considerations of the previous section will be our starting point to study nonholonomic integrators, 
and our equations will be conceptually equivalent to the proposed for systems with external forces 
(see mi). In the particular case of mechanical systems with linear constraint in the velocities, we 
study a subclass of our family of nonholonomic integrators with the property of preservation of the 
original nonholonomic constraints. 


2 Geometrical formulation of nonholonomic systems 


Let Q be a n-dimensional differentiable manifold, with local coordinates (q 1 )- The tangent bundle 
TQ, with induced coordinates (■ q l ,q l ), is equipped with two fundamental geometrical objects 1331: 
the Liouville vector field A and the vertical endomorphism S. In natural bundle coordinates we 
have 


A = q l 


d_ 

dq l 


S = dq l 


d_ 

dq l 


Consider a Lagrangian system, with Lagrangian L : TQ —» M, subject to nonholonomic constraints, 
defined by a submanifold D of the velocity phase space TQ. We will assume that dimH = 2 n — m 
and that D is locally described by the vanishing of m independent functions <j> a (the “constraint 
functions”). 

In geometrical terms, D’Alembert’s principle (or Chetaev’s principle for nonlinear constraints) 
implies that the constraint forces, regarded as 1-forms on TQ along D , take their values in the 
subbundle S*(TD° ) of T*TQ, where TD° denotes the annihilator of TD in T*TQ. In an intrinsic 
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way, the equations of motion can be written as (see E3H3) 


(ixv L -dE L )\ D G S*(TD°), 

X\ D G TD, 

where lol is the Poincare-Cartan 2-form defined by ujl = — d(S*(dL )) and El = A (L) — L is the 
energy function. 

In the sequel we will also assume that the following admissibility condition holds 


dimTD° = dimS*(TD°). 

This essentially means that the matrix (dcj) a /dq l ) has rank m everywhere. 

We now turn to the Hamiltonian description of the nonholonomic system on the cotangent bundle 
T*Q of Q 0 -21., 333- The canonical coordinates on T*Q are denoted by ( q l ,pt ), and the cotangent 
bundle projection will be ttq : T*Q —> Q. Assuming the regularity of the Lagrangian, we have that 
the Lagrangian and Hamiltonian formulations are locally equivalent. If we suppose, in addition, that 
the Lagrangian L is hyperregular, then the Legendre transformation Leg : TQ —> T*Q , (q l ,q l ) i—> 
( q l ,Pi = dL/dq l ), is a global diffeomorphism. The constraint functions on T*Q become = 
<f> a o Leg —1 , i.e. 

* a (q\Pi) = nq\ o-), 

OPi 


where the Hamiltonian H : T*Q 
(q , -7—), then 

OPi 


is dehned by H = El o Leg 1 . Since locally Leg 1 (q\pi ) = 


H = p t q l - L(q\ q v ) , 

where g* is expressed in terms of q l and pi using Leg~ 1 . 

The equations of motion for the nonholonomic system on T*Q can now be written as follows 

t)H 


Q = 
Pi = 


dpi 

dH m °\j 

dq l a dpj ^ ’ 


(5) 


together with the constraint equations d> a (q,p) = 0, where Tiij are the components of the inverse 
of the matrix (W J ) = (d 2 H / dpidpj). Note that 

,dd! a d 4> a 

(^« Jt )(M) = (^-°h9 )(q,p). 


The symplectic 2-form u>l is related, via the Legendre map, with the canonical symplectic form loq 
on T*Q. Let M denote the image of the constraint submanifold D under the Legendre transfor¬ 
mation, and let F be the distribution on T*Q along M, whose annihilator is given by 

F° = Leg*(S*(TD°)). 
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Observe that F° is locally generated by the m independent 1-forms 


u a = —— Hijdq 3 , 1 < a <m. 

dpi 

The “Hamilton equations” for the nonholonomic system can be then rewritten in intrinsic form as 


(ixUQ-dH)\ M G F° 

X\m £ TM. 


( 6 ) 


Suppose in addition that the following compatibility condition F 1 - P\TM = {0} holds, where “ T ” 
denotes the symplectic orthogonal with respect to ujq. Observe that, locally, this condition means 
that the matrix 

. /<9T a d'S> b \ 

ic) = (*"«<*) (7) 
is regular. On the Lagrangian side, the compatibility condition is locally written as 


det(C ab ) = det ) #0, 


( 8 ) 


where W %:> are the entries of the Hessian matrix 


d 2 L 


The compatibility condition is 


dq'd&S l<ij<n 

not too restrictive, since, taking into account the admissibility assumption, it is trivially verified by 
the usual systems of mechanical type (i.e. with a Lagrangian of the form kinetic minus potential 
energy), where the Ttij represent the components of a positive definite Riemannian metric. The 
compatibility condition guarantees in particular the existence of a unique solution of the constrained 
equations of motion © which, henceforth, will be denoted by Xh,m on the Hamiltonian side and 
£l,d on the Lagrangian side. 

Moreover, if we denote by Xjj the Hamiltonian vector field of H , i.e., Ix h ojq = dH then, using the 
constraint functions, we may explicitely determine the Lagrange multipliers A a as 


A q = -C ab X H (^ b 


Next, writing the 1-form 

. d^! a 

A = -C ab X H (y b )—- Hjidq 1 , 
d pj 

the nonholonomic equations are equivalently rewritten as 


q l = 


dH 


dpi 

dH 4 

* = ~W~ ‘ 


(9) 


for initial conditions (qo,po) G M and A = A idq 1 . We also denote by A = Leg*(A) the 1-form on 
TQ wich represents the constraint force once the Lagrange multipliers have been determined. 


Now, consider the flow F b : M —► M, t G / C K of the vector field Xh,m, solution of the 
nonholonomic problem. 









Since 0 is geometrically rewritten as 


= dH + A , 

(i% L d lol = dEi + A, with A = Leg*A., on the Lagrangian side) then 

l x hm 0q = d(ix HtM 0Q - H) - A , 


or, equivalently, 

Lx - H ,M d Q = d ( L ° Le 5 _1 ) - A • 

Now, from the dynamical definition of the Lie derivative, we have 

and integrating, we obtain the following expression, with some abuse of notation, 


— On = d 


f 


L o Ft dt) — 


F t * A 


( 10 ) 


where Ft is the flow of the vector field £ l,d • In next sections, we will study geometric integrators 
which verify a discrete version of equation m- 


3 “Generating functions” and nonholonomic mechanics 

Next, we will follow similar arguments for the construction of generating functions for symplectic 
or canonical maps |I . However, because of equation ©, we have that the nonholonomic flow is 
not a canonical transformation; i.e., 

FZu, Q -u,Q = d(J h Ft a) . (11) 

This description will allow us to construct a new family of nonholonomic integrators for equations 
©• Denote by nt : T*Q x T*Q —> T*Q, i = 1,2, the canonical projections. Consider the following 
forms 


0 = n* 2 e Q - 7rie Q , 

D = TT^COQ - TT^OJq = -d@ . 


Denote by iF h ■ Graph (T),) <—> T*Q x T*Q the inclusion map and observe that Graph(T)j) C M x M. 
Then, from m 


* F h ^ “ ( 7 r l|Graph(F h )) (F^loq w Q ) 

rh 


= frl I Graphs))* 


d 


F t * A 
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or, from m, 


*Fh 0 = (’ri|Graph(F h ))* 


fh 


L o Ft dt 



Let (qo,Po,qi,Pi) be coordinates in T*Q x T*Q in a neighborhood of some point in Graph(F/j). If 
{qo,Po->quPi) £ Graph(F/j) then ^“(gojPo) = 0 and 'S> a (qi,pi) = 0. Moreover, along Graph(T)j), 
qi = qi{qo,Po) and Pi = pi(qo,po), 


Pi dqi - p 0 dq 0 


d[ / L(q(t),q(t))dt 



Mq(t),q(t)), 


( 12 ) 


where (q(t),q(t)) = F t (q 0 ,q 0 ) with Leg(q 0 ,qo) = (qo,Po)- Here, F t denotes the flow of £ L>D . Equa¬ 
tion m is satisfied along Graph(i'X). 

Assume that, in a neighborhood of some point x G Graph (Fh), we can change this system of 
coordinates to a new coordinates (qo ? Qi ) - Denote by 

S h (qo,qi)=[ L(q(t),q(t))dt, 

J 0 

where q[t) is a solution curve of the nonholonomic problem with g(0) = q and q(h) = q\. This 
solution always exists for adequate values of go and q\ . In fact, observe that 

OH 

qi = Qo + h—(q 0 ,Po ) + o(h 2 ) , 

hence, since det (^ Jp f y p ^j / 0 , we locally have that po = po(qo, qi, h). But, in addition, (qo,Po) £ M; 
therefore ip a (qo,q\,h) = '& a (qo,Po(qo,qi,h)) = 0 . Then, the curve 


(q(t),q(t)) = Leg 1 (F t (q 0 ,po(qo,qi,h))) , 


verifies the required assumptions if <p a (qo,qi, h) = 0. 

Thus, we deduce that 1 

dS h f h ~ dq 

Po = — + / A (q(t), q(t))-^~ , 

oqo Jo vqo 

< 

dS h f h ~ <9g 

where (qo,qi) verifies the constraint functions <p a (qo,qi,h) = 0, now explicitely defined by 

dS h f h ~ da 

p a (qo,qi,h) = ^ a (qo,-- 7 ^(qo,qi) + A (q(t),q(t))-^-) , 1 <a<m, 


(13) 


(14) 


with q(t) solution of the nonholonomic problem with g(0) = go and q(h) = qy- 

x For a function f(x,y) with x,y £ ffi" we use the notation df/dx (respectively, df/dy) to write the partial 
derivative with respect the first n-variables (resp., the second n- variables). 
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Next, we will show how the group composite law of the flow F k 


FNh = F h ° ■ ■ ■ ° Fh 

N 

is expressed in terms of the corresponding “generating functions” S h . Moreover, the following 
Theorem will result in a new construction of numerical integrators for nonholonomic mechanics 
when we change the “generating function” and the constraint forces by appropriate approximations. 
As a generalization of Theorem 1 1.1 1 we have the following 

Theorem 3.1 The function S Nh , the “generating function” for Fj^h, is given by 

N -1 

S Nh (qo,qN) = sh (Qk,qk+ 1 ) , 

k =0 

where qk, 1 < k < N — 1, are points verifying 

D 2 S h (q k -x,q k ) + DiS h {q k ,qk+i) = A (q(t),q(t))-^-+ J A (q(t),q(t))-^j- , (15) 

and q(t) is a solution curve of the nonholonomic problem with q(0) = q k -i and q{h) = q k (respec¬ 
tively, q(h) = q k and q(2h) = q k +i) for the first integral (resp., second integral) of the right-hand 
side. 

Proof: It is suffices to prove the result for N = 2; that is, 

S 2h (qo,q 2 ) = S h (q 0 ,qi) +S h (qi,q 2 ) , 
where q\ verifies condition m- 
Since 

pi dqi - p 0 dq 0 = 

P 2 dq 2 ~ pi dqi = 

then 

P 2 dq 2 -Podqo = d (s h (q 0 ,qi) + S h (qi,q 2 )'j - J A (q(t), q{t)) - J^ A (q(t),q(t)) . 

Since the variables qi do not appear on the left-hand side term, it follows that 

0 = D 2 Si(qo, qi) + DiS!f(qi, q 2 ) — f A (q(t),q(t))-p-- [ A (q(t),q(t))-p- , (16) 

Jo dqi J h dq 0 

and for a choice of qi verifying m then 

S 2h (qo, qi) = S h (q 0 , qi) + S h (qi,q 2 ) 


dS h (q 0 ,qi)~ A (q(t),q(t)) 

Jo 


r'2h 


dS (qi,q 2 )~ / A(q(t),q(t)) 

Jh 




f2 h 
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is a “generating function of the first kind” of Foh because 


P 2 dq 2 - Po dq 0 



A (q(t),q(t)) . | 


Equations m determine a implicit system of difference equations which permit us to obtain q 2 from 
the initial data 90 and q \. An interesting consequence is that these equations preserve the constraint 
submanifold determined by the constraints p a = 0, 1 < a < m. In fact, if p a (qo, qi,h) = 0 (that is 
^ a (qo,Po ) = 0) then 

f)Q h fh „ F) n 

<p a (qi,q2,h) = * a ( qi ^(qo,qi)- / A(q(t), q(t))P~) , 

oqi Jo oqi 

and now applying m we obtain that 

P a (qi,q2,h) = ^“(<21,Pi) = 0 , 
since Fh(qo,po) = (qi,pi) and the flow preserves the constraints. 

The next remark will be a key result for the construction of nonholonomic integrators. 


Remark 3.2 Replace equation (EH) by 


dS h 

po = + a oO?o,gi), 

dS h ° h 

Pi = -w-«1 (?o, 9i) , 

oqi 


(17) 


where S h is a function of (qo, 9 i) coordinates and a h = cJq dqo + a\ dq\ and replace the constraints 
functions by 

dS h 

<P a (qo,qi,h) = ^ a (qoi~-Q^o + a o(9o,9i)) , (18) 

that is, 

Pi dqi — po dqo = dS h — a h , 


along (p a = 0 . 
Assume that 


det 


/ d 2 S h 


dan 


/o 


(19) 


y dqodqx dq x 

then, applying the implicit function theorem we have that, locally, q\ = 91 ( 90 , Po), and then the 
mapping 

G h (qo,Po) = (qi,Pi) 

is well-defined. 

Consider the mapping G^h defined by 

GNh = Gh o • • • 0 Gh ■ 


N 


12 










Following a similar argument to Theorem 13.II Grapl^G/v/J is described by 

dS Nh 


Po = 


PN = 


-(<?o,<?iv) + aQ h {qo,qN) 


dqo 

dS Nh 

dqN 


(qo, qN ) - ai h (q 0 , qjy) , 


( 20 ) 


where S Nh (q 0 ,q N ) = Ylk=o S h {q k ,q k +i) and a Nh (q 0 ,q N ) = Ylk=o ah (qk,qk+i)- Here, the q k % 
1 < k < N — 1, verify 

D 2 S h (q k -i,q k ) + D 1 S h (q k ,q k+ i) = ot\{q k _i,q k ) + 9fc+i), 1 < k < N - 1 . (21) 


3.1 Constraint error analysis 

As we have seen, if our “generating function” is S h , then we have exact preservation of the con¬ 
straints ip a . We now investigate what happens when the “generating function” is an approximation. 
We follow similar arguments to those in subsection 2.3.1 in m- 

Assume that Q , and also TQ and T*Q, are finite-dimensional vector spaces with inner product 
(.,.) and corresponding norm j| ||. 

Consider an “approximated generating function” S h and an approximated discrete constraint force 
a h = a!? dq l for the nonholonomic problem both of order r; hence, there exists an open set U C D 
with compact closure and constants c,di> 0, 1 < i < n, and H > 0 such that 

S h (qo,qi) = S h (q 0 ,qi)+C{q 0 ,q 1 ,h)h r+1 


oJj = 


fh _ 

/ A i{q(t), q(t)) dt + A(<?o, qi, h)h r+1 
Jo 


( 22 ) 

(23) 


for all solution q(t) of the nonholonomic problem with q( 0) = qo, q{h) = q\ and initial condition 
belonging to U and h < H. Here C and D{, 1 < i < n, are functions such that ||C((7o> qi, h) || < c 
and \\Di(q 0 ,qi,h)\\ < di on U. 

Taking derivatives we have that 

dS h 


a®, ( ®' ?l) = ■g^ {qo ’ qi) + wJ q °' quhW+1 


and also 


a, 


o(<?o,<?i) = (ao)iTT- = [ Mq{t),q(t))^-dt+ ^2^-(q 0 ,q ll h)h r+1 
dqo Jo oq 0 dq 0 


where now a h = a (j dqo + a\ dq\ 
Therefore, we deduce that 


dS 


<P a (qo,qi,h) = ^ a (qo,~g^Q +a 0 (qo,qi)) 

= * a (qo,—Q^ + J Q Hq(t),m)-^- o )+E a (q 0 ,qi,h)h r+1 
= ^ a (qo,Po) + E a (q 0 ,q 1 ,h)h r+1 = E a (q 0 ,q 1 ,h)h r+1 
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where E a are bounded functions. Then, the discrete algorithm preserves the constraints up to 
order r. 


3.2 Local error analysis 


Assuming that 


det 


d 2 S h 

dqodqi 


dqi ) 


7 ^ 0 , 


(24) 


we obtain a discrete flow G h : V C M —> M. It is easy to show, from conditions HU and HU, 
that Gh is an integrator of Xh,m of order r, following similar arguments to those used in the 
subsection above (see also Theorem 2.3.1., in mr 


4 Nonholonomic integrators 

In the sequel and for simplicity assume that Q is a vector space. Since S h (qo, q±) = jJ 1 L(q(t),q(t)) dt , 
where q(t) is a nonholonomic solution with q( 0) = qo and q(h) = q\ , using R.emark 13.21 we can 
obtain nonholonomic integrators by taking adequate approximations of the “generating function” 
S h and the extra-term A(q(t), q(t)). 

Consider, for instance, the approximation 

Sa(qo, qi) = hL(( 1 - a)q 0 + aqi, Ql ^ g ° ) , (25) 

for some parameter a £ [0,1]. (In general, we will write S^(qo, q\) ss S h (qo, q\ )-) 

A natural approximation of the constraint forces adapted to our choice of approximation for S h 
are 

J A (q(t),q(t))-^- « (1 - a)hA((l - a)q 0 + aqi, Ql ^ g ° ) , 

J A (q(t),q(t))-^- ~ ah\((l - a)q 0 + aqi, Ql ^ 90 ) . 


Consequently, equations m give us the following numerical method for nonholonomic systems 


D 2 Sa(qk-i,qk) + D 1 Sa{qk,qk+i) = ahA((l - a)q k -i + aq k , Qk - ) 
+(1 - a)hA((l - a)q k + aq k+ 1 , ——) , 1 < fc < IV — 1 , 


with 


initial condition satisfying 


^ a (<?o, qi,h) = ^“(<70, 9i) + (1 - a)&A((l - «)^o + a?i, gl h g ° )) = 0 
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Remark 4.1 Obviously, it is possible to produce a wider variety of discrete methods. For example, 


1 oh i 1 oh 


on _ _ on I _ on 

‘- , sym,a '2° a ' a ’ 

gives a second-order method for any a € [0,1]. Also, higher-order approximations of the function 
S h may be considered. 


Example 4.2 Nonholonomic particle. 

Consider the Lagrangian L : TR 3 —>■ R 

1 


L = -(x 2 + y 2 + z 2 )-(x 2 + y 2 ), 


subject to the constraint 


(f> = z — yx = 0 . 

It is easy to compute the nonholonomic differential equations 

2x + yxy 


x = — - 


1 + y 2 
y = -2 y 

—2 xy + xy 


z = 


1 + 2 / 


2 ’ 


where now the constraint 1-form is 


X 2 xy-xy 

A = ——— y-(dz - ydx) . 
l+y 2, 


Taking 


h 


S 1/2 (xo,y 0 ,z 0 ,x 1 ,yi,z 1 ) = - 


xi xq \ 2 f yi-yo\ 2 , fzi-zo 


h 


+ 


h ) 


+ 


h 


Xo + Xi \ 2 fyo + 2/1 N 2 


we obtain the nonholonomic integrator 


X\ - Xo , X\ + Xq X 2 X\ X 2 + X\ 

— h --- h- 


h 2 h " 2 

(ai+soXi/i+j/o) _ (xi—xo)(yi—yo) 

2 P 


l + (ai+M) 2 


.. , .. (x 2 +Xl)(j/2+yi) (x 2 -Xi)(j/2-yi) 

2/1 + 2/0 + 2-P- ! 2/2 + 2/1 


1+ (^+w) 2 


2/1 - 2/0 , 2/1 + 2/0 2/2 - 2/1 , 2/2 + 2/1 n 

--- h --- h -= 0 , 

h 2 h 2 

21 - z 0 z 2 Zl 


h h 

' (xi+a:o)(i/i + i/o) _ (a;i-a:o)(i/i-yo) (x2+x 1 )(y 2 +yi) _ (x 2 -x 1 )(y 2 -yi) 

2 P , 2 P 


h 

2 


1 + (yi+M) 2 


l+(^i) 2 
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The constraint function on M 3 x M 3 is 


t P a (x 0 ,yo,z 0l x 1 ,y 1 ,z 1 ,h) = 


, (xi+x 0 )(yi+yo) (a:i — a: 0 ) (?/l—J/0) 

Zl - ~0 _ » 2 ~ h? 

h 2 1 + (Ei+lfi>) 2 


+ 


yo 


- x 0 
h 


+ h 


Xi + 
2 


1 (x 1 +x 0 )(yi+yo) (xi—xo)(yi—yo) .. , _ 

h 2 ~ ti> . 2/1 + yo 

2 i + (ai+M) 2 2 


The following two hgures show the preservation of energy as a key point of comparison of compu¬ 
tational implementations of the method exposed above to other methods. 

The first figure compares the method introduced here to the traditional Runge-Kutta method of 
fourth order, showing an improvement in several orders of magnitude. Observe that, in this scale, 
the value of the energy in each step of our algorithm is practically undistinguishable from the initial 
value of the energy. 


energy 



The second figure is a comparison between our method and the one appeared in mm- A similar 
behaviour is observed. Nevertheless, a slightly better behaviour can also be appreciated, where the 
proposed algorithm shows on average a better preservation of the original energy. 
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energy 



For the same initial conditions and data, the following graph shows a very good behaviour of the 
constraint function evolution with time (notice the small scale). 


0 50 100 150 200 250 300 350 400 450 500 


5 Mechanical systems with linear constraints. Geometric numer¬ 
ical methods preserving constraints 


Suppose that the mechanical system, given by the Lagrangian L : TQ 


L (v q ) = ^g{v q ,v q ) ~V(q) 
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is subjected to nonholonomic constraints cj) a : TQ —> R, 1 < a < m. Since the nonholonomic 
constraints usually found in mechanics are linear in the velocities we will assume that 

4> a {qA) = 1 <a<m. 

From a geometric point of view, these linear constraints are determined by prescribing a distribution 
V on Q of dimension n — m such that the annihilator of V is locally given by 

V° = (//“ = nfdq 1 ; 1 < a < m). 

In this manner, the solutions of the nonholonomic Lagrangian system satisfy 


Vc(t)c(i) = -grad V(c(t)) + A (c(t)), c(t) G V c{t) , (26) 

where A is a section of V 1 - along c, and V ^ stands for the orthogonal complement of V with respect 
to the metric g. 

Since g is a Riemannian metric, the m x m matrix ( C ab ) = ( gfg ij ) is symmetric and regular. 
Therefore, we can explicitly determine 

aWW,«'(*)) = c ab z b (27) 

where ( C a b ) is the inverse matrix of ( C ab ) and the vector field Z a is defined by 


g(Z a ,Y) = g a {Y), for all vector field Y, 1 < a <m , 

that is, Z a is the gradient of the 1-form fi a . Thus, V 1 - = (Z a ), 1 < a < m. In local coordinates, 
we have 


d 

r/a _ jij , ,a 

z ~ s “-W 

By using the metric g and the distribution V we can obtain two complementary projectors 


V : TQ 
Q : TQ 


V, 


with respect to g. The projector Q is locally described by 


Q = C ab Z a ® / 


Using these projectors we can obtain the equations of motion as follows. A curve c(t) is a motion 
for the non-holonomic system if it satisfies the constraints, say, (j) a (c(t )) = 0, for all a, and, in 
addition, the “projected equation of motion” 

c(t )) = —U(grad V(c(t))) (28) 

is fulfilled. But these conditions are equivalent to 

c(t) € V c(t) , Vc( t )c(i) = —T’(grad V(c(t))) , 
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where V is the modified linear connection defined by 


Vjh = V X Y + (V A -Q)(y) 


for all vector fields X and Y on Q. 

Since the constraints are linear then, from m 

f) Qh fh _ f) n 

-Vi(qo)g tJ (qo)—r{qo,qi) + tf(qo)g l:, (qo) A(q(t),q(t))—- = o , i<a<m, 

oq 0 d % 

or, in terms of projectors, 

Q\g 0 ( y DiS h (q 0 ,qi)'j) = Q\ qo (di A(q(t),q(t))^j 

Moreover, the dynamics preserves the constraints 'h a which implies that 


( 29 ) 


(30) 


(qi, ^r—(qo,qi) - [ A(g(f),g(i))-|^) = 0 , 
oqi Jo oq\ 

or, in other words, 

Q\ qi (D 2 S h (q 0 , gi )) = Q\ qi (d 2 j A (q(t),q(t)j (31) 

Therefore, equations ob and OB show that the preservation of the exact constraints is equivalent 
to give a prescription about the relationship between the “generating function” and the constraint 
forces. 


Thus, equations m 

D 2 S h (q k -i,q k ) + D 1 S h (q k ,q k+ i) = J A(q(t),q{t))-^-+ A(q{t),q(t))-^- , 

can be rewritten using expression OB as follows 

V \<ik [piS h {q k - i,g fc )) + DiS h (q k ,q k+1 ) = V\ qk {^ j A(q(t),q(t))-^j + A (q(t),q(t))-^- , 

Q1 h g °(32) 

Now, considering an approximated generating function S h and an approximate constraint force 
a h = «o( 90 ) Q\ ) dqo + a ((go, qi) dqi, as in Remark 13.21 from the previous discussion, we now substi¬ 
tute the approximated constraint force by: 


a h = aQ(qo,qi)dq 0 

+'P\qMi(qo,qi) dqi) + Q\ qi (p 2 S h {qo,qi)^j) 


Therefore for S h and a h equations OB are rewritten as 

V \qk { D 2 S h (q k -i,q k )j + Di§ h (q k , q k+1 ) = V\ qk (a^ (q k -i, q k )^j + <4(g fc , q k + 1 ), (33) 


19 




for 1 < k < N — 1. The importance of equations (ET31) is that they generate an algorithm which 
automaticaly preserves the exact constraint functions <h a . In fact, if we apply the projector Q to 
Equations (PI) we obtain: 

Q\q k (DiS h (q k ,q k + 1 )) = Q\ qk ( ao(q k ,q k +i ) (34) 

or 

floh 

<P a (qk, Qk+i, h) = ^ a (q k , 

dqo 

that is, the constraints are satisfied. 

Therefore the geometric algorithm that we have obtained work as follows: 

V\ qk (D2S h {q k -i,qk)) + DiS h {q k ,q k+1 ) = V\ qk (oq (g fc _i, %)) + a%(q k , q k+1 ), 

with initial condition satisfying: 

<p a (qo,qi,h) = 0 

Choosing ckq and a\ in T>°, we obtain equations for nonholonomic integrators with more geometric 
flavour: 

Geometric nonholonomic integrator 

V\ qk ^D 2 S h (q k -i,qk) + D 1 S h {q k , q k +i)^ = 0 

which is interpreted as a discretization of Equations Pt 

Vc(t)c(i) = —'P(grad (V(c(i))) 

In a future work we will study from numerical and geometrical points of view this particular subclass 
of geometric integrators. 


(qk,Qk+i) +Uo(q k ,qk+i)) = 0 


5.1 Nonholonomic integrators preserving constraints 

For the class of integrators introduced in Section 0J we find the following family of nonholonomic 
integrators preserving constraints: 

V \q k faSaiqk-i, %)) + D i s a(qk, qk+ 1 ) = ahV\ qk ^A((l - a)q k -i + aq k , Qk ^ k ~ l )\ 

+(1 - a)/iA((l - a)q k + aq k+ 1, gfc+1 ^ gfc ) ; 1 < k < N - 1 , 

with initial condition satisfying 

~tf(qo)9 ij (qo)^(qo, ?i) + (1 - a) h^(q 0 )g^(q 0 )A J ((l - a)q 0 + aq u = 0 . 

d% h 
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Example 5.1 (The nonholonomic particle revisited) 


1 / Xi - X 0 _ ^ Xi + x 0 \ _ x 2 - X 1 _ h x 2 + Xi 2/1 { Z\ — Z() 


1 + 2/f V h 


h 

2 


1_ (xi+x 0 )(y 1 +y 0 ) _ (x 1 -x 0 )(y 1 -y 0 ) 

2 h ? 2/1 + 2/0 


1+^T ' 


i + y\ \ h 

(■ x2+xi)(y 2 +yi) _ (x2-xi)(y2-yi) 


1 + (M+ffi) 2 


+ 


~W 


2/2 + 2/1 


1+ (^) 2 


(x 1 +x 0 )(y 1 +y 0 ) (xi-x 0 )(y 1 -y 0 ) 

Vl 2 ~ 7? 

yi+yol 2 


1 + 2/2 i + (3O+S0) 


2/1 - 2/0 , 2/1 + 2/0 2/2 - 2/1 , 2/2 + 2/1 n 

- ; - h --- h -= 0 , 

h 2 h 2 

y\ ( z\ - z 0 \ z 2 Z\ 2/1 ^ ^1 ~ _ ^i + ®o 


1 + 2/? V h 

h 


h 1 + 2/f V ^ 


2 (xi+xp)(yi+yo) _ (x 1 -x 0 ){y 1 -y 0 ) (x 2 +x 1 )(y 2 +yi) _ (x 2 -xi)(y 2 -yi) 

Vl 2 h 2 _|_2 ft 3 


1 + 21? i+(ai+M) 2 

(a:i+xo)(yi+3/o) (.xi—xo)(yi—yo) 


1+ (^+W) 2 


2/1 


T 3- 


2/1 +2/o 


1 + 2/2 l + (ai+m) 2 


with initial condition satisfying 


r (xi+ao)(yi+yo) (xi-x 0 )(yi-yo) 

~a/ j\ 2-1 -S'O 2 h? 

* = - 1 t - 

xi - io , t n + xo h (£i±a+±a) _ fc£fa»! 9l + 90 

+ “ ~ + ft —2“ “2- 1 + (E +a+ -2“ 


For the same initial conditions and data, the following graph shows the exact preservation of the 
constraint function evolution with time of our algorithm. 


Constraint preservation 


i r~ 


i i i 
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6 Conclusion 


A new numerical algorithm has been proposed for nonholonomic mechanics. This algorithm is 
based in the underlying geometry of nonholonomic systems. For mechanical systems with linear 
constraints, a geometric integrator preserving constraints is proposed. 

In future work, we will explore reduction schemes for discrete systems using the approach of gen¬ 
erating functions. It is also interesting to use generating functions of different kinds; in a recent 
work EH, we have shown that generating functions of second class generate algorithms which are 
symplectic (in some sense) for discrete optimal control theory (see also EH). Moreover, we may 
easily extend the generating function technique in order to consider variable time stepping and also 
the time-dependent case and it would be possible to use this formalism for classical field theories. 
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